Carga de librerías

library(pacman)
packages = c("MASS","knitr","tidyverse","car",'dplyr','kableExtra',"tidyr","readr","magrittr","VIM","GGally","igraph","ggplot2","sf","leaflet","webshot2","ggmosaic","corrplot","readxl","stringr","ggpubr")
pacman::p_load(char=packages)

Carga de ficheros

Para crear el dataset con el que vamos a tratar en este proyecto, hemos extraido varios archivos de la web del portal de datos abiertos del Ayuntamiento de Valencia. En ellos tenemos diferente información acerca de los 88 barrios que hay en Valencia, como pueden ser el número de zonas verdes, precio del alquiler, actividad comercial, renta, etc. Antes de atacar las preguntas que nuestro conjunto resolverá, vamos a cargar los datos y unirlos en un único dataset, con una variable común para todos, el barrio.

Vulnerabilidad

El primer dataset Vulnerabilidad nos da información general del barrio, como la densidad de población, la renta media, o el estado de vulnerabilidad. Esta última variable será de gran interés en nuestro análisis posterior.

vuln <- read_delim("./data/vulnerabilidad.csv", 
    delim = ";", escape_double = FALSE, col_types = cols(`Geo Point` = col_skip(), 
        `Geo Shape` = col_skip(), `Densitat_p` = col_skip()), 
    trim_ws = TRUE)%>%arrange(nombre)

colnames(vuln)[colnames(vuln) == "nombre"] <-"Barrio"
colnames(vuln)[colnames(vuln) == "Index_Gl_1"] <-"Indice_Vuln"


vuln$Barrio<-factor(vuln$Barrio,levels = unique(vuln$Barrio))
vuln$Indice_Vuln<-factor(vuln$Indice_Vuln,levels = c("Vulnerable","Pot. Vulnerable","No Vulnerable"))

head(vuln)

Podemos ver en el código que es importante transformar la variable Barrio en un factor, para poder graficar y tratar la información de forma adecuada. Repetiremos este proceso en cada conjunto de datos, además de poner a todos el mismo nombre para poder unirlos más adelante.

Precio de compra y alquiler

Estos dos datasets de compra y alquiler nos presentan informaciones similares, que es la media de precios de compra y alquiler en nuestros barrios en los años 2022 y 2010. Ya que son conjuntos muy similares, vamos a adelantarnos al próximo paso y fusionarlos en un único dataset, llamado precios.

p_compra <- read_delim("./data/precio-de-compra-en-idealista.csv", 
    delim = ";", escape_double = FALSE, col_types = cols(`Geo Point` = col_skip(), 
        `Geo Shape` = col_skip(), Fecha_creacion = col_skip(),`Max_historico (Euros/m2)` = col_skip(), 
        Año_Max_Hist = col_skip()), 
    trim_ws = TRUE)%>%arrange(BARRIO)

p_compra$BARRIO<-factor(p_compra$BARRIO,levels = unique(p_compra$BARRIO))

p_alquiler <- read_delim("./data/precio-alquiler-vivienda.csv", 
    delim = ";", escape_double = FALSE, 
    trim_ws = TRUE, col_types = cols(`Geo Point` = col_skip(), `Geo Shape` = col_skip(), Fecha_creacion = col_skip(),`Max_historico (Euros/m2)` = col_skip(), 
        Año_Max_Hist = col_skip()))%>%arrange(BARRIO)

p_alquiler$BARRIO<-factor(p_alquiler$BARRIO,levels = unique(p_alquiler$BARRIO))

precios<-full_join(p_compra,p_alquiler,by="BARRIO",suffix = c(" de compra"," de alquiler"))
colnames(precios)[colnames(precios) == "BARRIO"] <-"Barrio"

head(precios)

Recibos IBI

Vamos ahora con el dataset IBI, que nos da información de los diferentes recibos del IBI (Impuesto sobre Bienes Inmuebles) entre los años 2021 y 2023. Este conjunto nos va a dar una muy buena visión acerca de la actividad del barrio, tanto comercial como cultural, turística, religiosa, industrial, etc.

Debido a que en ningún momento vamos a tratar con tiempo en este dataset, vamos a eliminar los años haciendo la media de las observaciones de cada barrio durante estos tres años, para así obtener tantas observaciones como barrios, ya que si no habrá conflictos a la hora de unir los datos.

ibi <- read_delim("./data/rebuts-ibi-2022.csv", delim = ";", escape_double = FALSE, col_types = cols(ID = col_skip(), geo_shape = col_skip(), geo_point_2d = col_skip(), Año = col_skip()), trim_ws = TRUE)%>%
  arrange(Barrio)%>%
  mutate_at(vars(-all_of(c("Distrito","Barrio"))), ~as.numeric(sub(",",".",.)))

ibi$Barrio<-factor(ibi$Barrio,levels = unique(ibi$Barrio))

# Hacemos la media de las observaciones de cada barrio en los tres año y nos quitamos 2/3 de las observaciones
ibi <- ibi %>% group_by(Barrio) %>%  mutate(across(where(is.numeric), mean, na.rm=TRUE))%>%distinct()

head(ibi)

Bancos por barrio

Por último, vamos con nuestro último conjunto de datos, barrios, que contiene mucha información acerca de la ubicación de las entidades bancarias en nuestra ciudad. Debido a que nosotros solo vamos a tratar con barrios y no con direcciones ni nada similar, hemos decidido que lo más interesante de este conjunto es el número de bancos que podemos encontrar en cada barrio (puede ser un buen indicador de riqueza o pobreza). Guardaremos esta información en un nuevo dataset llamado num_bancos.

bancos <- read_delim("./data/bancs-en-via-publica-bancos-en-via-publica.csv", 
    delim = ";", escape_double = FALSE, col_types = cols(gid = col_skip(), 
        `Num. Policia` = col_skip(), geo_point_2d = col_skip()), 
    trim_ws = TRUE)%>%arrange(Barrio)

bancos$Barrio%<>%gsub("[0-9] - ","",.)
bancos$Barrio%<>%factor(levels = unique(bancos$Barrio))

num_bancos<-bancos%>%group_by(Barrio)%>%summarize(Num_bancos=n())
head(num_bancos)

Población

Este conjunto contiene el area y la población de los diferentes barrios. Hemos considerado este dataset debido a que los datos de areas y densidad de población que nos proporcionaba el conjunto vulnerabilidad y barrios no cuadraban con nuestras búsquedas y carecían de sentido. Por ello hemos considerado este otro que se ajusta mucho mejor. Para calcular la densidad usaremos la función mutate.

Este dataset se trata de un directorio donde se encuentran un total de 176 archivos Excel. Dos de estos aportan información relacionada con Valencia, mientras que el resto dan información específica de cada barrio (por cada barrio hay dos archivos, uno en castellano y otro en valenciano). Por tanto, hemos creado un programa que abra uno a uno los archivos en un idioma y obtenga la información que deseamos para después guardarla en un dataframe.

Para ello, primero hemos abierto uno de los archivos Excel (ver Figura 1) y se ha estudiado como abordar el problema. Después, se ha cargado uno de ellos, y se ha buscado la manera de obtener la información referida al nombre del barrio, el área y la población del mismo (recuadros amarillos en la Figura 1). Figura 1: ejemplo de uno de los archivos Excel analizados.

#   SE HA PUESTO CON eval = FALSE PORQUE TARDA MUCHO EN CORRER
# Obtener las rutas de todos los archivos con los que se va a trabajar
archivos_barrios <- list.files(path = "./data/Barrios2022", pattern = 'Districte', full.names = TRUE)

# Creamos las listas donde vamos a guardar los datos de interés (la longitud es de 87 porque hay información sobre 87 barrios)
area <- 1:87
poblacion <- 1:87
nombre <- 1:87

# Un bucle que recorre todos los archivos que utilizaremos y extrae los datos de interés
for (i in 1:length(archivos_barrios)) {
  Barrio <- read_excel(archivos_barrios[i], sheet = "Padrón")
  fila <- grep('Padró' ,Barrio[[1]])[1]
  nombre_barrio <- Barrio[[1]][fila]
  nombre[i] <-trimws(toupper(substring(nombre_barrio, 36)))
  fila <- grep('Superficie' ,Barrio$...2)
  area[i] <- Barrio$...2[fila + 1]
  fila <- grep('Personas', Barrio[[1]])[1]
  poblacion[i] <- Barrio[[1]][fila + 1]
}

demografico <- data.frame(nombre, area, poblacion)

# Hay algunos nombres que se han guardado mal, por lo que vamos a tener que tratarlos. Vemos que todos los que estan mal comienzan por ".  " o por "[0-9].  "
for (i in 1:length(demografico$nombre)) {
  if (grepl('^[^A-Za-z]',demografico$nombre)[i]) {
    demografico$nombre[i] <- str_remove(demografico$nombre[i], "([0-9]|). ")
  }
}

# También falta el dato de un barrio, asi que hay una fila de más
demografico <- demografico[!(demografico$nombre == '88'), ]

# Vemos quales de los nombres del dataframe nuevo no están en el principal
precios$Barrio[!(demografico$nombre %in% precios$Barrio)]

# Hay que poner de la misma maera que estan en el data-frame pricipal
demografico$nombre[demografico$nombre == "GRAN VIA"] <- "LA GRAN VIA"
demografico$nombre[demografico$nombre == "EXPOSICIÓ"] <- "EXPOSICIO"
demografico$nombre[demografico$nombre == "CIUTAT UNIVERSITÀRIA"] <- "CIUTAT UNIVERSITARIA"
demografico$nombre[demografico$nombre == "SANT MARCEL·LÍ"] <- "SANT MARCEL.LI"
demografico$nombre[demografico$nombre == "CAMÍ REAL"] <- "CAMI REAL"
demografico$nombre[demografico$nombre == "MONT-OLIVET"] <- "MONTOLIVET"
demografico$nombre[demografico$nombre == "FONTETA DE SANT LLUÍS"] <- "LA FONTETA S.LLUIS"
demografico$nombre[demografico$nombre == "CIUTAT DE LES ARTS I DE LES CIÈNCIES"] <- "CIUTAT DE LES ARTS I DE LES CIENCIES"
demografico$nombre[demografico$nombre == "EL CABANYAL- EL CANYAMELAR"] <- "CABANYAL-CANYAMELAR"
demografico$nombre[demografico$nombre == "EL BOTÀNIC"] <- "EL BOTANIC"
demografico$nombre[demografico$nombre == "BETERÓ"] <- "BETERO"
demografico$nombre[demografico$nombre == "CAMÍ FONDO"] <- "CAMI FONDO"
demografico$nombre[demografico$nombre == "CIUTAT JARDÍ"] <- "CIUTAT JARDI"
demografico$nombre[demografico$nombre == "LA BEGA BAIXA"] <- "LA VEGA BAIXA"
demografico$nombre[demografico$nombre == "CAMÍ DE VERA"] <- "CAMI DE VERA"
demografico$nombre[demografico$nombre == "ORRIOLS"] <- "ELS ORRIOLS"
demografico$nombre[demografico$nombre == "SANT LLORENÇ"] <- "SANT LLORENS"
demografico$nombre[demografico$nombre == "CASES DE BÀRCENA"] <- "LES CASES DE BARCENA"
demografico$nombre[demografico$nombre == "MAUELLA"] <- "MAHUELLA-TAULADELLA"
demografico$nombre[demografico$nombre == "\t\nBORBOTO"] <- "BORBOTO"
demografico$nombre[demografico$nombre == "\t\nBENIMAMET"] <- "BENIMAMET"
demografico$nombre[demografico$nombre == "EL CASTELLAR-L'OLIVERAR"] <- "CASTELLAR-L'OLIVERAL"

# Guardar el RData
save(demografico, file = "./data/demografico.RData")

Una vez creado, vamos a cargarlo:

load("./data/Demografico.RData")

colnames(demografico)[colnames(demografico) == "nombre"] <-"Barrio"

demografico$Barrio<-factor(demografico$Barrio,levels = unique(demografico$Barrio))

demografico%<>%
  arrange(`Barrio`)%>%
  mutate(across(-c("Barrio"), as.numeric))%>%
  mutate(`Densidad`=`poblacion`/`area`)

head(demografico)

Hemos tenido que aplicar una transformación a las columnas a traves de un mutate y un across ya que estas eran de tipo carácter, y no de tipo numérico.

Fusion de los dataset

df<-vuln%>%full_join(demografico,by="Barrio")%>%full_join(num_bancos,by="Barrio")%>%full_join(precios,by="Barrio")%>%full_join(ibi,by="Barrio")

dim(df)
[1] 94 63
tail(df)

Vemos como a la hora de fusionar todos los datos en un solo dataset, tenemos un problema, y es que contamos con más observaciones de las esperadas. Deberíamos tener un total de 88 observaciones (una por cada barrio), pero en cambio, tenemos 92. Mirando el final del dataset, vemos como efectivamente hay cuatro observaciones que no se corresponden con lo deseado, así que vamos a arreglarlo.

# El barrio MONTOLIVET se llama únicamente en el dataset "vuln" y "num_bancos" como MONT-OLIVET
levels(vuln$Barrio)[vuln$Barrio=="MONT-OLIVET"]<-"MONTOLIVET"
levels(num_bancos$Barrio)[num_bancos$Barrio=="MONT-OLIVET"]<-"MONTOLIVET"

# Lo mismo ocurre con la Fonteta de sant lluis y el dataset "ibi
levels(ibi$Barrio)[ibi$Barrio=="FONTETA DE SANT LLUIS"]<-"LA FONTETA S.LLUIS"

# Además, la primera y última observación de barrios no corresponden a ningún barrio, por lo que tenemos que eliminarlas
num_bancos%<>%slice(-c(1,length(num_bancos$Num_bancos)))

Vemos como dos de estas incongruencias se debían a la distinta forma de escribir el nombre de los barrios, mientras que las otras dos simplemente se debían a que algún dataset contenía información de barrios desconocidos, lo cual es mejor eliminar directamente.

Una vez solucionado, volvemos a crear el dataset:

df<-vuln%>%full_join(demografico,by="Barrio")%>%full_join(num_bancos,by="Barrio")%>%full_join(precios,by="Barrio")%>%full_join(ibi,by="Barrio")

dim(df)
[1] 90 63
head(df)

Con esto, ya tenemos el número de observaciones deseadas, asi que podemos proceder con el siguiente paso.

Selección de variables inicial

Observando las 69 variables con las que cuenta nuestro conjunto, vemos como claramente hay muchas que no necesitamos. Primero, tenemos todos los códigos de los barrios, que prácticamente cada dataset de los anteriores contaba con una o más variables de este estilo, y con distintos nombres entre sí. Vamos a empezar eliminandolas aplicando expresiones regulares, ya que todas cuentan con una cosa en común, que empiezan por “cod”:

codigos<-grepl("^[Cc]od",colnames(df))
df%<>%select(-colnames(df)[codigos])

Con esto nos hemos quitado un total de 11 variables, pero aún podemos hacer más. También tenemos otra variable redundante, que son los distritos. Como el distrito de compra es el único que no tiene ningún valor perdido, usaremos ese, y además, lo transformaremos en factor:

distritos<-colnames(df)[grepl("^DISTRITO|Distrito",colnames(df))]
distritos<-distritos[distritos!="DISTRITO de compra"]

df%<>%select(-distritos)

colnames(df)[colnames(df) == "DISTRITO de compra"] <-"Distrito"

df$Distrito<-factor(df$Distrito,levels = unique(df$Distrito))

df%<>%relocate(Barrio,Distrito,Indice_Vuln)

Por último, vemos como dentro del dataset IBI tenemos por un lado las variables que indica el número de recibos de un cierto tipo y en otra el importe. Consideramos que nos pueden ser más de utilidad las segundas, y para no ser reduntantes vamos a eliminar las de número de recibos. Además, algunas de estas cuentan con entradas decimales, lo cual es un poco extraño para lo que la variable representa.

num<-grepl("^Num\\.",colnames(df))
df%<>%select(-colnames(df)[num])

Finalmente tenemos nuestro conjunto de datos cargado y liberado de variables innecesarias, vamos a echar un vistazo:

dim(df)
[1] 90 36
head(df)
numeric_cols <- sapply(df, is.numeric)
factor_cols <- sapply(df, is.factor)

# Crear un data frame para el gráfico de barras
summary_df <- data.frame(
  variable = colnames(df),
  tipo = ifelse(numeric_cols,"Numeric", "Factor"),
  count = c(sum(numeric_cols), sum(factor_cols))
)

# Crear un gráfico de barras
ggplot(summary_df, aes(x = tipo, y = count, fill = tipo)) +
  geom_bar(stat = "identity") +
  labs(title = "Número de columnas numéricas y factoriales",
       x = "Tipo de variable",
       y = "Número de columnas") +
  theme_minimal()

ggplot(summary_df, aes(x = variable, y = tipo, fill = tipo)) +
  geom_bar(stat = "identity") +
  labs(title = "Número de columnas numéricas y factoriales",
       x = "Tipo de variable",
       y = "Número de columnas") +
  theme_minimal()+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1,size=rel(0.8)))

Tenemos un total de 33 variables numéricas (cuantitativas) y 3 variables de tipo factor (cualitativas). Una vez creado y depurado nuestro conjunto de forma preliminar, vamos a plantear las preguntas que queremos resolver y acabar de poner nuestro dataset a punto.

Preguntas a resolver

Claramente este conjunto de datos gira entorno a un objeto, los barrios de Valencia. Además, la variable del Índice de Vulnerabilidad consideramos que tiene muchísimo interés. Por tanto, todo nuestro trabajo va a consistir en encontrar las variables que más influyen en esta vulnerabilidad, y encontrar la manera en la que estas contribuyen.

Hemos tratado de responder estas preguntas, primero, realizando un estudio de la correlación entre todas las variables numericas de nuestro dataset, así como de las dos variables categoricas del mismo, creando así unos grupos dentro de las variables. Posteriormente, hemos realizado diferentes gráficas para poder relacionar la compleja relación entre las mismas.

Estudio de la correlación para variables numericas

df_numeric <- select_if(df, is.numeric)

# Correlación de Pearson
cor_pearson <- cor(df_numeric, method = "pearson", use = "complete.obs")

# Correlación de Spearman
cor_spearman <- cor(df_numeric, method = "spearman", use = "complete.obs")
# Pearson
cor_pearson_long <- as.data.frame(cor_pearson) %>%
  rownames_to_column(var = "Variable1") %>%
  gather(key = "Variable2", value = "Correlation", -Variable1)

# Filtrar las correlaciones mayores a 0.8 y menores de 1 (para positivos y negativos)
strong_correlations <- cor_pearson_long %>%
  filter(abs(Correlation) > 0.8, abs(Correlation) < 1) %>%
  filter(!duplicated(t(apply(.[, c("Variable1", "Variable2")], 1, sort))))

# Spearman
cor_spearman_long <- as.data.frame(cor_spearman) %>%
  rownames_to_column(var = "Variable1") %>%
  gather(key = "Variable2", value = "Correlation", -Variable1)

# Filtrar las correlaciones mayores a 0.8 y menores de 1 (para positivos y negativos)
strong_correlations_spearman <- cor_spearman_long %>%
  filter(abs(Correlation) > 0.8, abs(Correlation) < 1) %>%
  filter(!duplicated(t(apply(.[, c("Variable1", "Variable2")], 1, sort))))

print(strong_correlations)
                                  Variable1                                Variable2 Correlation
1                                 poblacion                               atur_16_64   0.8877622
2                                Index_Soci                               renda_mitj   0.9089642
3          Precio_2022 (Euros/m2) de compra                               renda_mitj   0.9132029
4                                Index_Glob                               Index_Soci   0.9131050
5          Precio_2022 (Euros/m2) de compra                               Index_Soci   0.8495091
6  Imp.Recibos Actv.Almacen-Estacionamiento           Importe Recibos personalidad F   0.9171733
7            Imp. Recibos Actv. Residencial           Importe Recibos personalidad F   0.9974711
8                   Importe Recibos totales           Importe Recibos personalidad F   0.8846799
9              Imp. Recibos Actv. Comercial           Importe Recibos personalidad J   0.8890860
10              Imp. Recibos Actv. Oficinas           Importe Recibos personalidad J   0.8068381
11                  Importe Recibos totales           Importe Recibos personalidad J   0.8616861
12              Imp. Recibos Actv. Oficinas         Importe Recibos sin personalidad   0.8422134
13           Imp. Recibos Actv. Residencial Imp.Recibos Actv.Almacen-Estacionamiento   0.9223414
14                  Importe Recibos totales Imp.Recibos Actv.Almacen-Estacionamiento   0.9276058
15                  Importe Recibos totales             Imp. Recibos Actv. Comercial   0.8503477
16          Imp. Recibos Actv. Espectaculos              Imp. Recibos Actv. Cultural   0.8718010
17                  Importe Recibos totales           Imp. Recibos Actv. Residencial   0.8869359
print(strong_correlations_spearman)
                                  Variable1                                Variable2 Correlation
1                                risc_pobre                               renda_mitj  -0.8511416
2                                Index_Soci                               renda_mitj   0.9697542
3                                Index_Glob                               renda_mitj   0.8989763
4          Precio_2022 (Euros/m2) de compra                               renda_mitj   0.8801309
5                                Index_Soci                               risc_pobre  -0.8745461
6                                Index_Glob                               risc_pobre  -0.8313238
7                                Index_Glob                               Index_Soci   0.9495787
8          Precio_2022 (Euros/m2) de compra                               Index_Soci   0.8411217
9            Importe Recibos personalidad F                                poblacion   0.8010323
10 Imp.Recibos Actv.Almacen-Estacionamiento           Importe Recibos personalidad F   0.9362086
11             Imp. Recibos Actv. Comercial           Importe Recibos personalidad F   0.9045753
12           Imp. Recibos Actv. Residencial           Importe Recibos personalidad F   0.9954510
13                  Importe Recibos totales           Importe Recibos personalidad F   0.9530750
14             Imp. Recibos Actv. Comercial           Importe Recibos personalidad J   0.8096754
15              Imp. Recibos Actv. Oficinas           Importe Recibos personalidad J   0.8337853
16                  Importe Recibos totales           Importe Recibos personalidad J   0.8357099
17             Imp. Recibos Actv. Comercial Imp.Recibos Actv.Almacen-Estacionamiento   0.9526551
18           Imp. Recibos Actv. Residencial Imp.Recibos Actv.Almacen-Estacionamiento   0.9356137
19                  Importe Recibos totales Imp.Recibos Actv.Almacen-Estacionamiento   0.9526551
20           Imp. Recibos Actv. Residencial             Imp. Recibos Actv. Comercial   0.9041554
21                  Importe Recibos totales             Imp. Recibos Actv. Comercial   0.9482460
22                  Importe Recibos totales           Imp. Recibos Actv. Residencial   0.9562593
corrplot(mar = c(1, 0, 1, 0), cor_spearman, method = "circle", type = "upper", order = "hclust",
         tl.col = "black", tl.srt = 45,
         tl.cex = 0.3,
         cl.lim = c(-1, 1), cl.cex = 0.3, cl.ratio = 0.1,
         addCoef.col = "black", number.cex = 0.25, coef.cex = 0.3,
         col = colorRampPalette(c("#6D9EC1", "white", "#E46726"))(200),
         main = "Correlaciones Spearman")

En este grafico de correlaciones de Spearman podemos observar visualmente como de correlacionadas están todas nuestras variables entre todas ellas, ya que estamos representando las correlaciones de spearman, éstas no tienen por qué ser relaciones lineales. Ya que no se puede observar de forma clara las relaciones, realizaremos un diagrama de Red donde poder ver solo las más grandes tanto de Pearson como de Spearman.

Primero, calcularemos las correlaciones de spearman que no estén en pearson, ya que una correlación alta de pearson implica tambien una correlación alta de spearman, pero saber que es de pearson nos da una información mayor a si sabemos que es de spearman, ya que nos dice que la relación es lineal. Posteriormente, representaremos en un diagrama de red las correlaciones de pearson y las correlaciones de spearman que no sean también de pearson.

# Unir ambos data frames
all_correlations <- rbind(strong_correlations, strong_correlations_spearman)

# Crear una función para ordenar las variables (para que las parejas sean consistentes en orden)
order_pair <- function(df) {
  df %>%
    rowwise() %>%
    mutate(Variable1 = pmin(Variable1, Variable2),
           Variable2 = pmax(Variable1, Variable2)) %>%
    ungroup()
}

# Ordenar pares en ambos data frames
all_correlations_order <- order_pair(all_correlations)
strong_correlations_order <- order_pair(strong_correlations)

# Restar ambos dataframes
spearman_not_in_pearson <- anti_join(all_correlations_order, strong_correlations_order, 
                                 by = c("Variable1", "Variable2"))
# Asegurarse de que no hay enlaces de un nado a si mismo
spearman_not_in_pearson <- spearman_not_in_pearson %>%
                            filter(Variable1 != Variable2)
set.seed(137)

spearman_not_in_pearson <- spearman_not_in_pearson %>%
  dplyr::select(Variable1, Variable2, Correlation) %>%
  mutate(Method = "Spearman")

strong_correlations <- strong_correlations %>%
  dplyr::select(Variable1, Variable2, Correlation) %>%
  mutate(Method = "Pearson")

# Combina los datos de Pearson y Spearman
combined_correlations <- rbind(strong_correlations, spearman_not_in_pearson)

# Crea un grafo a partir de los datos combinados
graph_data <- graph_from_data_frame(combined_correlations, directed = FALSE)

# Ajustar atributos del nodo
V(graph_data)$color <- "skyblue"
V(graph_data)$size <- 6
V(graph_data)$frame.color <- "black"

# Ajustar atributos de las aristas
cor_min <- 0.8
cor_max <- 1.0
width_min <- 1
width_max <- 5
E(graph_data)$width <- ifelse(E(graph_data)$Method == "Pearson",
                              (abs(E(graph_data)$Correlation) - cor_min) / (cor_max - cor_min) * (width_max - width_min) + width_min,2)

E(graph_data)$lty <- ifelse(E(graph_data)$Method == "Pearson", 1, 3)

layout <- layout_with_fr(graph_data)

par(mar = c(0, 0, 1.5, 1))
plot(graph_data, layout = layout, vertex.label.color = "black", vertex.label.cex = 0.7,
     vertex.label.dist = 1.2,
     edge.label = NA,
     edge.color = ifelse(E(graph_data)$Correlation < 0, "red", "gray"),
     main = "Red de Correlaciones")

legend("topright",
       legend = c("Pearson > 0.8", "Pearson < -0.8", "Spearman > 0.8", "Spearman < -0.8"),
       col = c("gray", "red", "gray", "red"),
       lty = c(1,1,3,3),
       cex = 0.8)

En este diagrama queda muy claro las relaciones entre variables con correlaciones altas, diferenciando correlacion pearson de la de spearman, y correlaciones negativas de positivas. Podemos observar que los importes de los distintos recibos están muy correlacionados y estos a su vez lo están con la población a través de los recibos de personalidad física, la población además ésta correlacionada con el desempleo. Por otro lado, los recibos de actividades de espectaculo y de actividades de cultura están muy correlacionados. Por último, las variables que reflejan el nivel adquisitivo también están muy correlacionadas entre ellas e inversamente correlacionadas con el riesgo de pobreza. De aquí podríamos si quisieramos reducir mucho el número de variables de nuestro dataset, elegir una variable de cada grupo mostrado y eliminar el resto, pues aportan poca información.

Estudio de la correlación para variables categoricas

Realizamos un test chi-cuadrado para conocer si hay una relacion clara entre el distrito en el que se encuentra el barrio y el indice de vulnerabilidad de ese barrio. El resultado es un p-value de 0.01966, por lo que podemos rechazar la hipotesis nula de que las variables son independientes a un nivel de significancia del 0.05 (nivel de significancia estándar).

chisq_result_distrito <- chisq.test(df$Distrito, df$Indice_Vuln)
print(chisq_result_distrito)

    Pearson's Chi-squared test

data:  df$Distrito and df$Indice_Vuln
X-squared = 55.568, df = 36, p-value = 0.01966

Podemos graficar en un mosaico estas dos variables para observar que efectivamente existe una relación entre ambas.

library(ggmosaic)
ggplot(data = df) +
    geom_mosaic(aes(weight = 1, x = product(Distrito), fill = Indice_Vuln)) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, size = rel(0.8))) +
    labs(title = "Mosaic Plot de Distrito vs Indice_Vuln", x = "Distrito", y = "Índice de Vulnerabilidad")+ scale_fill_manual(values = c("#F8766D", "#619CFF", "#00BA38", "#7F7F7F"))

Detección de anomalías

Detección de valores perdidos

Antes de tratar con nuestros datos, vamos a analizar la situación de nuestro dataset. Primero, vamos a analizar los valores perdidos.

aggr(df, prop = FALSE, combined = TRUE, numbers = TRUE, sortCombs = TRUE)

Este gráfico nos muestra las observaciones con valores perididos y en qué columnas se hayan. Por ejemplo, para la primera observación, vemos como todas las columnas a excepción de dos cuentan con un NA, y así sucesivamente, hasta llegar a ver que hay 70 observaciones sin ningún valor perdido.

sum(is.na(df))
[1] 210

Podemos ver que la cantidad de valores perdidos en nuestro conjunto no es precisamente pequeña, y principalmente se debe al hecho de que no todos los conjuntos de datos que hemos fusionado contenían información de todos los barrios, por lo que a la hora de unirlos todos se han generado NAs en las observaciones donde no existían datos.

Una cosa que salta a la vista de las variables es esa observación que cuenta con casi todos los valores perdidos, que es la del barrio “RAFALELL-VISTABELLA”, que no cuenta con ninguna información numérica en nuestro dataset. Por ello, lo mejor que podemos hacer es eliminarla.

df[df$Barrio=="RAFALELL-VISTABELLA",]
df%<>%drop_na("Importe Recibos personalidad F")

colSums(is.na(df))
                                  Barrio                                 Distrito 
                                       0                                        0 
                             Indice_Vuln                               Zones verd 
                                       4                                        4 
                              turismes_e                               atur_16_64 
                                       4                                        4 
                              renda_mitj                               risc_pobre 
                                       4                                        4 
                              Index_Equi                               Index_Soci 
                                       4                                        4 
                              Index_Glob                                     area 
                                       4                                        2 
                               poblacion                                 Densidad 
                                       2                                        2 
                              Num_bancos         Precio_2022 (Euros/m2) de compra 
                                       2                                       17 
        Precio_2010 (Euros/m2) de compra       Precio_2022 (Euros/m2) de alquiler 
                                      17                                       17 
      Precio_2010 (Euros/m2) de alquiler           Importe Recibos personalidad F 
                                      17                                        0 
          Importe Recibos personalidad J         Importe Recibos sin personalidad 
                                       0                                        0 
Imp.Recibos Actv.Almacen-Estacionamiento             Imp. Recibos Actv. Comercial 
                                       0                                        0 
             Imp. Recibos Actv. Cultural             Imp. Recibos Actv. Deportiva 
                                       0                                        0 
      Imp.Recibos Actv.Edificio singular          Imp. Recibos Actv. Espectaculos 
                                       0                                        0 
           Imp. Recibos Actv. Industrial      Imp.Recibos Actv.Obras urbanizacion 
                                       0                                        0 
      Imp.Recibos Actv.Ocio y Hosteleria              Imp. Recibos Actv. Oficinas 
                                       0                                        0 
           Imp. Recibos Actv. Religiosas           Imp. Recibos Actv. Residencial 
                                       0                                        0 
Imp.Recibos Actv.Sanidad y Beneficiencia                  Importe Recibos totales 
                                       0                                        0 

En esta tabla podemos ver que variables son las que cuentan con datos perdidos, y por tanto las que debemos procesar. Vemos como las variables del dataset Vulnerabilidad cuentan todas con 4 NAs, que se deben a las observaciones de los barrios que no cuentan con Indice de Vulnerabilidad. Como esta va a ser nuestra variable objetivo durante todo el trabajo, hemos preferido no tratar con ellos, para no alterar así el resultado de nuestras observaciones.

Además, dado que algunas columnas, como las de precios, tienen un gran número de NAs, podemos usar el estudio de correlaciones que hemos visto anteriormente para sustituir el valor perdido por el equivalente en una de las columnas correlacionadas (usando una regresión, por ejemplo). En caso de tener un NA en la columna correlacionada, usaremos la mediana de los casos de su misma vulnerabilidad. En el caso de los precios de compra de 2022, usaremos la renta media del barrio, que cuentan con una correlación de 0.91:

reg<-lm(`Precio_2022 (Euros/m2) de compra`~renda_mitj,df)
  
df%<>%mutate(`Precio_2022 (Euros/m2) de compra`=ifelse(is.na(`Precio_2022 (Euros/m2) de compra`)&!is.na(renda_mitj),reg$coefficients[1]+renda_mitj*reg$coefficients[2],`Precio_2022 (Euros/m2) de compra`))

df %<>%
  group_by(`Indice_Vuln`) %>% 
  mutate(`Precio_2022 (Euros/m2) de compra`=ifelse(is.na(`Precio_2022 (Euros/m2) de compra`),median(`Precio_2022 (Euros/m2) de compra`,na.rm = TRUE),`Precio_2022 (Euros/m2) de compra`))%>%
  ungroup() 

Dado que el resto de variables con NAs no cuentan con ningun índice de correlación extremadamente altos, imputaremos los datos faltantes con la mediana de su mismo grupo de vulnerabilidad.

df %<>%
  group_by(`Indice_Vuln`) %>% 
  mutate(`Precio_2010 (Euros/m2) de compra`=ifelse(is.na(`Precio_2010 (Euros/m2) de compra`),median(`Precio_2010 (Euros/m2) de compra`,na.rm = TRUE),`Precio_2010 (Euros/m2) de compra`))%>%
  ungroup() 
df %<>%
  group_by(`Indice_Vuln`) %>% 
  mutate(`Precio_2022 (Euros/m2) de alquiler`=ifelse(is.na(`Precio_2022 (Euros/m2) de alquiler`),median(`Precio_2022 (Euros/m2) de alquiler`,na.rm = TRUE),`Precio_2022 (Euros/m2) de alquiler`))%>%
  ungroup() 
df %<>%
  group_by(`Indice_Vuln`) %>% 
  mutate(`Precio_2010 (Euros/m2) de alquiler`=ifelse(is.na(`Precio_2010 (Euros/m2) de alquiler`),median(`Precio_2010 (Euros/m2) de alquiler`,na.rm = TRUE),`Precio_2010 (Euros/m2) de alquiler`))%>%
  ungroup() 
df %<>%
  group_by(`Indice_Vuln`) %>% 
  mutate(`Num_bancos`=ifelse(is.na(`Num_bancos`),median(`Num_bancos`,na.rm = TRUE),`Num_bancos`))%>%
  ungroup() 

Veamos las conclusiones:

aggr(df, prop = FALSE, combined = TRUE, numbers = TRUE, sortCombs = TRUE)

Tras esto hemos logrado pasar de tener un número muy elevado de NAs a tener solo 4 en ciertas varibles. Estos NAs están en los barrios que carecen de índice de vulverabilidad, por lo que tendríamos que esperar a ponerles una etiqueta a estos barrios para librarnos de los NAs de forma adecuada.

Detección de outliers

Vamos a tratar los outliers de nuestro conjunto antes de empezar a trabajar.

Para la detección de outliers vamos a usar los métodos 3-sigma y boxplot, con las funciones definidas en la práctica 5.

reglasigma<-function(x){
  x<-x[!is.na(x)& is.numeric(x)]
  out <- logical(length(x)) 
  
  for(i in 1:length(x)){
    if(abs(x[i]-mean(x))>3*sd(x)){
      out[i]<-TRUE
    }
  }
  if (all(!out)){
    return(NA)
  } else {
    return(out)
  }
}

reglaboxplot<-function(x){
  x<-x[!is.na(x)& is.numeric(x)]
  out <- logical(length(x)) 
  
  for(i in 1:length(x)){
    if(x[i]>quantile(x,0.75)+1.5*IQR(x)){
      out[i]<-TRUE
    } else if (x[i]<quantile(x,0.25)-1.5*IQR(x)){
      out[i]<-TRUE
    }
  }
  if (all(!out)){
    return(NA)
  } else {
    return(out)
  }
}

Vamos a aplicar las funciones vistas en busca de posibles outliers:

outliers <- df %>%
  summarise(across(where(is.numeric), list(Sigma = ~sum(reglasigma(.)), Boxplot = ~sum(reglaboxplot(.))),.names = "{col};{fn}"))

outliers%<>%
  pivot_longer(cols=everything(), names_to = "Var",values_to = "Valor")%>%
  separate(Var,into=c("Variable", "Regla"),sep=";")%>%
  spread(key=Regla,value=Valor)

outliers

Viendo que la función boxplot detecta un número excesivo de outliers contando las pocas observaciones que tenemos, vamos a hacer caso a la regla sigma, y en caso de que haga falta modificar los outliers, solamente trataremos los que esta detecta, pasandolos a la mediana al igual que el ejemplo anterior, o usando alguna otra columna que esté muy correlacionada.

Observando y graficando detenidamente nuestro conjunto, hemos observado dos clases de outliers, un tipo del cual solo hemos detectado un caso, que concluimos que se debe a un fallo de inputación, y otro caso que se debe a un dato correcto, pero demasiado extremo, que rompe con la tendencia del resto de datos. Vamos a ver un ejemplo de cada uno.

Fallo de inputación

Este un ejemplo gráfico de otra forma de detectar ouliters. En el caso de la variable risc_pobre, el valor introducido para el barrio de Benimaclet, distaba 43.76 veces el rango intercuartílico de la mediana de la distribución.

Por tanto, se ha considerado un error de input y se le ha seleccionado un nuevo valor. Para ello, se ha tenido en cuenta que el análisis que se ha realizado ha sido mediante box-plots, donde se diferencian las distribuciones en función de la varaible categórica Indice_Vuln. Por tanto, para que no altere esta gráfica, el valor de la observación de Benimaclet se ha sustituido por la mediana correspondiente a la distribución con su misma vulnerabilidad.

# Corrección de outliers gráfica
ggplot(df, aes(x = Indice_Vuln, y = risc_pobre)) + geom_boxplot() + ggtitle('Distribución risc_pobre antes de tratar')

risc_pobre_filtrada <- df %>%
  filter(Indice_Vuln == 'Vulnerable') %>%
  filter(risc_pobre < 100) %>%
  select(risc_pobre)
df$risc_pobre[df['Barrio'] == 'BENIMACLET'] <- median(risc_pobre_filtrada[[1]])

ggplot(df, aes(x = Indice_Vuln, y = risc_pobre)) + geom_boxplot() + ggtitle('Distribución risc_pobre después de tratar')+ labs(x='Vulnerabilidad')

Vemos como era un outlier que no se correspondía con los datos, y al eliminarlo, tenemos una visión mucho más clara de nuestro conjunto.

Datos correctos

Un ejemplo de outlier puede verse en la variable que muestra la actividad cultural del barrio, viendo como la ciudad de las artes y las ciencias tiene un valor muchísimo más alto que el resto:

p<-ggplot(df[c("Imp. Recibos Actv. Cultural","Indice_Vuln")], aes(x = Indice_Vuln, y =.data[["Imp. Recibos Actv. Cultural"]])) +
    geom_boxplot() + 
    ggtitle(paste("Relación con Imp. Recibos Actv. Cultural")) + labs(x='Vulnerabilidad')
  print(p)

  print(df$`Imp. Recibos Actv. Cultural`[df$Barrio=="CIUTAT DE LES ARTS I DE LES CIENCIES"])
[1] 470265.7

Vemos como dentro de los barrios no vulnerables el de la ciudad de las artes y las ciencias tiene una actividad muchísimo mayor, con un valor de 470265.7. Aun así, debido a que este dato no se debe a un error a la hora de introducir el valor en el conjunto, pero se debe a que el barrio tiene una actividad cultural mayor debido a su situación. Por tanto, hemos decidido mantender estos outliers en nuestro dataset. En ciertos casos interesará eliminarlos siguiendo la regla sigma como hemos visto anteriormente, pero como nos limitamos a un análisis exploratorio, nos limitaremos a detectarlos.

#Analisis Univariante

Realizamos un estudio de las frecuencias en nuestras variables categoricas, que son dos variables de mucho interes en nuestro análisis, el indice de vulnerabilidad y los distritos.

frecuencias_distrito <- table(df$Distrito, useNA = "always")
df_frecuencias_distrito <- as.data.frame(frecuencias_distrito)
names(df_frecuencias_distrito) <- c("Distrito", "Frecuencia")

ggplot(df_frecuencias_distrito, aes(x = Distrito, y = Frecuencia)) +
  geom_bar(stat = "identity", fill = "blue") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=6), plot.title = element_text(hjust = 0.5)) +
  labs(title = "Frecuencia de Distritos", x = "Distrito", y = "Frecuencia")

frecuencias_Indice_Vuln <- table(df$Indice_Vuln, useNA = "always")
df_frecuencias_Indice_Vuln <- as.data.frame(frecuencias_Indice_Vuln)
names(df_frecuencias_Indice_Vuln) <- c("Indice_Vuln", "Frecuencia")

ggplot(df_frecuencias_Indice_Vuln, aes(x = Indice_Vuln, y = Frecuencia)) +
  geom_bar(stat = "identity", fill = "blue") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=8), plot.title = element_text(hjust = 0.5)) +
  labs(title = "Frecuencia de Indice Vulnerabilidad", y = "Frecuencia", x='Vulnerabilidad')

Vemos como la distribución no es uniforme entre todas las clases, sino que en el caso de los distritos cada uno tiene un número de barrios muy diferenciados, y en el caso de los indices de vulnerabilidad, más de la mitad de los barrios son no vulnerables mientras que el resto se dividen en vulnerables y pot. vulnerable. También observamos que no tenemos datos perdidos en los distritos y únicamente 4 valores perdidos en el índice de vulnerabilidad.

A continuación, estudiamos la variable densidad ya que es una medida de interés, podemos observar que a pesar de tener algun pico en alto como en densidades muy pequeñas, se distribuye bastante uniformemente. También calculamos algunos parametros que pueden ser de interés, como que la media de la densidad en Valencia es de 187.4078.

ggplot(df, aes(x = Densidad)) +
  geom_histogram(binwidth = 10,   
                 fill = "blue", 
                 color = "black") +
  labs(title = "Histograma de Densidad",
       x = "Densidad",
       y = "Frecuencia") +
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5))

summary(df$Densidad)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
  0.5426  76.0667 176.2565 187.4078 301.6097 529.1489        2 

#Analisis Bivariante

Al contrario de lo que podríamos pensar, en nuestra red de correlaciones no vemos una unión entre el precio de compra del metro cuadrado y el precio de alquiler del metro cuadrado, así que vamos a analizar esta relación. Vemos que la correlación entre el metro cuadrado de compra y de alquiler es alta pero no tanto como podríamos esperar.

cor_pearson_precios <- cor(df$`Precio_2022 (Euros/m2) de compra`, df$`Precio_2022 (Euros/m2) de alquiler`, method = "pearson", use = "complete.obs")
cor_pearson_precios
[1] 0.6443731
cor_spearman_precios <- cor(df$`Precio_2022 (Euros/m2) de compra`, df$`Precio_2022 (Euros/m2) de alquiler`, method = "spearman", use = "complete.obs")
cor_spearman_precios
[1] 0.651632

Esto implica que, en algunos barrios, el precio de compra puede ser relativamente alto mientras que el precio de alquiler puede ser bajo, o viceversa. Puede ser util si queremos buscar una vivienda para alquilar o comprar ya que vemos que barrios son mejores para cada caso.

Aprovecharemos este analisis para ver como afectan estas dos variables a nuestra variable de interés, la vulnerabilidad. Al añadirle como tercera variable el indice de vulnerabilidad de forma de color de punto, podemos observar que hay una relación entre el indice de vulnerabilidad y el precio del metro cuadrado de compra pero no así con el precio del metro cuadrado de alquiler. Esto quiere decir que los barrios donde la compra de vivienda es más cara son barrios no vulnerables, mientras que si el precio de alquiler es alto no querrá decir que el barrio no sea vulnerable.

ggplot(df, aes(x = `Precio_2022 (Euros/m2) de compra`, y = `Precio_2022 (Euros/m2) de alquiler`)) +
    geom_point() +
    labs(title = "Relación entre Precio por Metro Cuadrado de Compra y Alquiler",
         x = "Precio por Metro Cuadrado de Compra",
         y = "Precio por Metro Cuadrado de Alquiler")

ggplot(df, aes(x = `Precio_2022 (Euros/m2) de compra`, y = `Precio_2022 (Euros/m2) de alquiler`, color=Indice_Vuln)) +
    geom_point() +
    labs(title = "Relación entre Precio por Metro Cuadrado de Compra y Alquiler",
         x = "Precio por Metro Cuadrado de Compra",
         y = "Precio por Metro Cuadrado de Alquiler")

Estudio de la vulnerabilidad

En primer lugar, hemos representado un mapa de valencia, donde se han superpuesto formas geométricas (polígonos) que indican los distintos barrios de Valencia. Para ello, hemos cargado un archivo .geojson y hemos creado un nuevo datra-frame llamado datos_geojson. En este, además de guardar las variables relacionadas con el mapa interactivo, se han añadido las columnas que han resultado interesantes para representar.

#primero creamos un mapa con los barrios y distrito
library(sf)
library(ggplot2)

# Lee el archivo GeoJSON
datos_geojson <- st_read("./data/barris-barrios.geojson")
Reading layer `barris-barrios' from data source 
  `C:\Users\aitoo\Documents\AED\PROYECTO\ProyectoAED2023\ProyectoAED2023\data\barris-barrios.geojson' 
  using driver `GeoJSON'
Simple feature collection with 88 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -0.432535 ymin: 39.27893 xmax: -0.2753685 ymax: 39.56659
Geodetic CRS:  WGS 84
colnames(datos_geojson)[2] <- 'Barrio'

# Selecciono las columnas de df que quiero utilizar en la siguiente gráfica
columnas <- c('Barrio', 'Distrito', 'Indice_Vuln', 'area','poblacion','Densidad')

# Creo el dataframe que voy a utilizar para las gráficas
datos_geojson <- datos_geojson %>%
  full_join(df[columnas], by = 'Barrio')

Para crear los mapas, se ha utilizado la librería Leaflet, que permite la creación de mapas interactivos. En la siguiente gráfica, se muestran los barrios de Valencia, coloreados en función de la vulnerabilidad de los mismos. En primer luegar, nos ha parecido adecuado empezar representando cómo se distribuye la población de Valencia. Para ello, en un mapa, se han representado los distritos coloreados en función de la cantidad de personas que viven allí.

# vamos a hacer un mapa donde se muestre la cantidad de gente que vive en cada barrio
distritos_geojson <- st_read("./data/districtes-distritos.geojson")
Reading layer `districtes-distritos' from data source 
  `C:\Users\aitoo\Documents\AED\PROYECTO\ProyectoAED2023\ProyectoAED2023\data\districtes-distritos.geojson' 
  using driver `GeoJSON'
Simple feature collection with 22 features and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -0.432535 ymin: 39.27893 xmax: -0.2753685 ymax: 39.56659
Geodetic CRS:  WGS 84
# El nombre de uno de los distritos no coincide -> lo cambio a mano
distritos_geojson$nombre[distritos_geojson$nombre == 'POBLATS DE L\'OEST'] <- "POBLATS LOEST"
colnames(distritos_geojson)[2] <- 'Distrito'

# Selecciono las columnas que quiero utilizar en la siguiente gráfica
columnas <- c('Barrio', 'Distrito', 'Indice_Vuln', 'area','poblacion')

# Creo el dataframe que voy a utilizar para las gráficas
distritos_geojson <- distritos_geojson %>%
  full_join(df[columnas], by = 'Distrito')

# Creo una variable nueva que sea la poblacion del distrito donde se encuentra cada barrio
distritos_geojson <- distritos_geojson %>%
  group_by(Distrito) %>%
  mutate(Pob_dist = sum(poblacion, na.rm = T))
# Creo los popup del mapa
popups_dist <- paste0("<b>", distritos_geojson$Distrito, "</b>", "<hr>", paste('Población: \n', distritos_geojson$Pob_dist))

# Escojo una paleta de colores
pal <- colorBin("YlOrRd", domain = distritos_geojson$Pob_dist, bins = 6)

# Creo el mapa
leaflet(data = distritos_geojson) %>%
  addTiles() %>%
  addPolygons(fillColor = pal(distritos_geojson$Pob_dist),
              weight = 1,
              opacity = 1,
              highlightOptions = highlightOptions(color = "white",
                                                  weight = 2,
                                                  bringToFront = TRUE),
              color = 'black',
              fillOpacity = 0.8,
              popup = popups_dist) %>%
  addLegend(data = distritos_geojson,
            position = 'bottomright',
            pal = pal, values = ~Pob_dist,
            title = 'Población',
            opacity = 1)

Esta misma gráfica, la pondeomos hacer para los barrios de Valencia

# Creo los popup del mapa
popups <- paste0("<h3>", datos_geojson$Barrio, "</h3>", "<hr>", '<b>', 'Población: ', '</b>', datos_geojson$poblacion)

# Escojo una paleta de colores
pal <- colorBin("YlOrRd", domain = datos_geojson$poblacion, bins = 6)

# Creo el mapa
leaflet(data = datos_geojson) %>%
  addTiles() %>%
  addPolygons(fillColor = pal(datos_geojson$poblacion),
              weight = 1,
              opacity = 1,
              highlightOptions = highlightOptions(color = "white",
                                                  weight = 2,
                                                  bringToFront = TRUE),
              color = 'black',
              fillOpacity = 0.8,
              popup = popups) %>%
  addLegend(data = datos_geojson,
            position = 'bottomright',
            pal = pal, values = ~poblacion,
            title = 'Población',
            opacity = 1)

Se observa que, la mayoría de las personas viven en la zona céntrica, y que los barrios más poblados se encuentran en esta, a pesar de ser los barrios más pequeños. Esto lo podemos coroborrar también observándo la siguiente gráfica, donde se meuestra la densidad de población de cada barrio.

# Creo los popup del mapa
popups <- paste0("<h3>", datos_geojson$Barrio, "</h3>", "<hr>", '<b>', 'Población: ', '</b>', datos_geojson$Densidad)

# Escojo una paleta de colores
pal <- colorBin("YlOrRd", domain = datos_geojson$Densidad, bins = 6)

# Creo el mapa
leaflet(data = datos_geojson) %>%
  addTiles() %>%
  addPolygons(fillColor = pal(datos_geojson$Densidad),
              weight = 1,
              opacity = 1,
              highlightOptions = highlightOptions(color = "white",
                                                  weight = 2,
                                                  bringToFront = TRUE),
              color = 'black',
              fillOpacity = 0.8,
              popup = popups) %>%
  addLegend(data = datos_geojson,
            position = 'bottomright',
            pal = pal, values = ~Densidad,
            title = 'Densidad',
            opacity = 1)

Efectivamente, observamos que los barrios con una mayor densidad de población se encuentran en el centro, siendo la que mayor densidad tiene El Calvari, mientras que el barrio donde más gente vive es Benicalap.

En la siguiente figura, se han coloreado los barrios en función de su vulnerabilidad.

library(leaflet)
# Creo los popup del mapa
popups <- paste0("<h3>", datos_geojson$Barrio, "</h3>", "<hr>", datos_geojson$Indice_Vuln,'<br>', '<b>', 'Población: ', '</b>', datos_geojson$poblacion)

# Escojo una paleta de colores
pal <- colorFactor(c('red','gray','blue','green'), levels = levels(datos_geojson$Indice_Vuln))

# Creo el mapa
leaflet(data = datos_geojson) %>%
  addTiles() %>%
  addPolygons(fillColor = pal(datos_geojson$Indice_Vuln),
              weight = 1,
              opacity = 1,
              highlightOptions = highlightOptions(color = "white",
                                                  weight = 2,
                                                  bringToFront = TRUE),
              color = 'black',
              fillOpacity = 0.8,
              popup = popups) %>%
  addLegend(data = datos_geojson,
            position = 'bottomright',
            pal = pal, values = ~Indice_Vuln,
            title = 'Vulnerabilidad',
            opacity = 1)

Observando el mapa, vemos que los tanto El Calvari como Benicalap son considerados como barrios vulnerables. También se observa que en su mayoría, en las zonas colindantes a los barrios vulnerables, se encuentran más barrios vulnerables.

A continuación, vamos a seguir estudiando la vulnerabilidad de los barrios pero mediante diagramas de barras, el cual nos puede dar una idea más cuantitativa. En el siguiente diagrama, se muestran la cantidad de barrios que tiene cada distrito, coloreados en función de la vvulnerabilidad.

Sin embargo, para hacernos una idea de la cantidad de gente que vive en zonas vulnerables, es más adecuado el siguiente gráfico, donde se muestran la cantidad de personas que viven en cada distrito. Asimismo, al igual que en el anterior, se ha hecho una distinción en colores en función de si estas personas viven en barrios vulnerables.

Observamos claramente, que, la grafica anterior podía dar una idea errónea en cuanto a la cantidad de gente que vive en zonas vulnerables. Como muestra esta gráfica, la mayoría de personas en Valencia viven en barrios no Vulnerabes, aunque la cantaidad de personas que viven en barrios vulnerables o potencialmente vulnerables es mayor.

Dado que en este trabajo estudiamos la relación entre diferentes variables numéricas con la vulnerabilidad de un barrio, los Box-plot resultan la mejor opción para observar esta dependencia. En estos, somos capaces de estudiar la distribución de las variables numéricas en función de la vulnerabilidad. De esta manera, si las distribuciones que encontramos muy parecidas (ver siguiente figura), decimos que no se observa una relación clara (al menos cualitativamente), y por tanto consideraremos que esta variable no es influyente.

Por otro lado, existen algunas variables que sí han mostrado diferentes distribuciones en función de la vulnerabilidad. A continuación, se muestran las que nos han parecido de mayor interés. Las siguientes figuras muestran aquellas que nos han parecido esperables. Es decir, no es de extrañar que en los barrios vulnerables la tasa de paro sea mayor, o que el índice global (que mide la capacidad de un barrio en términos de innovación) sea mayor en los barrios no vulnerables. Asimismo, observamos que tanto la polbación como la densidad de población es mayor en barrios vulnerables. Por otro lado, en esta última, cabe destacar que la diferenciaentre la densidad de población en barrios vulnerables y la de barrios potencialmente vulnerables o no vulnerables.

Asimismo, nos parece interesante comparar los precios del m\(^2\) del año 2010 y 2022.

Vemos que para todas las distribuciones han sufrido un descenso de los precios, aunque para el caso de potencialmente vulnerable es el que menor cambio ha tenido.

Al estudiar todos los Box-plot realizado, ha habido otra gráfica que nos ha llamado la atención. Se trata de la variable Imp. Recibos Actv. Espectaculos (ver siguiente figura A), que parece que debido a outliers no se puede observar la distribución. Sin embargo, al representar el diagrama de barras, nos hemos dado cuenta de lo que realmente sucedía. A pesar de que en un principio pensáramos que podría tatarse de otro caso de error de input, al analizar estas gráficas concluimos que no era más que un outlier. Esta instancia de Imp. Recibos Actv. Espectaculos corresponde a la Ciudad de las Artes y Ciencias y, no es de extrañar que el importes de recivos relacionados con las actividades y espectáculos sea un outlier. Asimismo, vemos que para la mayoría de las instancias el valor de esta variable es de 0. Esto tampoco parece del todo realista, por lo que esta variable no parece del todo fiable y no se utilizará para realizar el análisis de la vulnerabilidad de barrios.

Por otro lado, al no tener una definición de lo que representa cada variable en la mayoría de datasets, hemos tenido que asumirlo a partir de los nombres. A pesar de ello, solo una no nos ha cuadrado con los datos esperados. Nosotros, hemos asumido que la variable turismes_e daba un valor proporcional al turismo de cada barrio. Así, cabría esperar que la ciudad de las Artes y las Ciencias sea la que mayor turismo tiene. Sin embargo, si vemos los barrios a los que corresponden los mayores valores de esta variable son:
Barrio turismes_e
CIUTAT JARDI 117.21
L’AMISTAT 76.11
L’ILLA PERDUDA 72.54
LA CARRASCA 22.85
LA VEGA BAIXA 45.20
NA ROVELLA 15.06

En un intento de contextualizar estos datos, hemos tratado de entender por qué estos barrios podrían ser los más turísticos, pero no hemos conseguido encontrar el por qué. Por tanto, hemos pensado que la carga del dataset se había hecho mal, por lo que hemos vuelto a descargar los datos y abierto el archivo, pero parece que los datos son correctos. Por tanto, concluimos que no entendemos del todo lo que representa esta variable.

Volviendo al estudio de la vulnerabilidad, en el siguiente gráfico de barras apliadas podemos ver la vulnerabilidad correspondiente al valor de atur_16_64 (relacionada con el índice de paro para personas entre 16 y 64 años). Como era de esperar, los valores más grandes correspondienden a los barrios vulnerables y a medida que disminuye el valor los barrios tienden a ser no vulnerables. El hecho de que se vea una clara segragación de la vulnerabilidad al ordenar el índice de paro nos indica que estas dos variables están correlacionadas.

Sin embargo, cuando estudiamos el índice global ocurre lo contrario. Los valores más grandes corresponden a barrios no vulnerables, y cuando nos moveos a valores inferiores los barrios vulnerables son mayoritarios. De nuevo, al representar las barras de manera ordenada, vemos zonas donde dominan los barrios vulnerables, zonas no vulnerables, y en medio las zonas potencialmente vulnerables. Por tanto, vemos que este índice también está relacionado con la vulnerabildiad de un barrio.

ggplot(df, aes(x = reorder(Barrio, -Index_Glob), y = Index_Glob, fill = Indice_Vuln)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("#F8766D", "#619CFF", "#00BA38", "#7F7F7F"), name = "Vulnerabilidad") +
  labs(title = "Barras apiladas") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(y="Index_Glob", x="Barrio") + 
  theme(axis.text.x = element_blank())

Para acabar, nos gustaría añadir que nos habría gustado haber realizado un análisis más avanzado con técnicas de cloustering para estudiar de una manera más amplia la correlación entre estas variables numéricas y la vulnerabilidad, pero es algo que se escapa de nuestro conocimiento. Por tanto, hemos intentado ver las variables que influyen en esta de manera cualitativa mediante las gráficas comentadas.

En conclusión, para este análisis hemos tratado de contextualizar las variables más relevantes, observarlas, y juzgar su bondad comparándolas con otros datasets o con nuestro conocimiento. Hemos conseguido estudiar la correlación entre las variables numéricas de una manera cuantitativa y sacar los grupos que mayor correlación tienen entre sí. Asimismo, mediante las gráficas de Box-plot y los diagramas concluimos que las variables que más afectan a este son: .

# scatter de risc_pobre - Index_soci
df%>%
      mutate(Indice_Vuln = ifelse(Indice_Vuln == 'Pot. Vulnerable', 'Vulnerable', Indice_Vuln)) %>%
      mutate(Indice_Vuln = ifelse(Indice_Vuln == '1', 'Vulnerable', Indice_Vuln)) %>%
      mutate(Indice_Vuln = ifelse(Indice_Vuln == '3', 'No Vulnerable', Indice_Vuln)) %>%
      mutate(Indice_Vuln = factor(Indice_Vuln, levels = c("Vulnerable", "No Vulnerable", "Pot. Vulnerable"))) %>%
ggplot(aes(x=risc_pobre, y = Index_Soci)) +
    geom_point(aes(color =Indice_Vuln)) + 
      ggtitle(paste("Relación entre risc_pobre y Index_Soci"))

Apéndice

Box-plot de tódas las variables numéricas en función de la vulnerabilidad de los barrios.